package CPStats;
use warnings;
use strict;
use Pod::Usage;
### The documentation of this module is at the __END__ of this file
## user modules
use Statistics;
use SimpleError;

my @comp_ops = (
    [ '>=', '<'  ], # mode = 2'b00
    [ '<=', '>'  ], # mode = 2'b01
    [ '>' , '<=' ], # mode = 2'b10
    [ '<' , '>=' ], # mode = 2'b11
);

my %discrete_desc = (
    'P'  => 'Pass',
    'F'  => 'Fail',
    '*F' => 'Fail',
    '*U' => 'Overflow Fail',
    '*L' => 'Underflow Fail',
);

my @comp_subs = map(
    [
        eval( sprintf('sub { $_[0] %s $_[1] }', $_->[0]) ),
        eval( sprintf('sub { $_[0] %s $_[1] }', $_->[1]) )
    ],
    @comp_ops
);


### class method
sub limit_desc {
    my ($ranges, $comp_mode, $na) = @_;
    my @desc;
    my $ops = $comp_ops[$comp_mode];
    if ( !( $comp_mode & 1) ) { # edge is on the left side
        push @desc, "$ops->[1] $ranges->[0]";
    }
    for (@$ranges) {
        push @desc, "$ops->[0] $_";
    }
    if ( $comp_mode & 1 ) { # edge is on the right side
        push @desc, "$ops->[1] $ranges->[ $#$ranges ]";
    }
    push @desc, 'n/a' if $na;
    return \@desc;
}

sub new {
    my $class = shift;
    my $self = +{ @_ };
    bless $self, $class;
}

sub set {
    my $self = shift;
    %{ $self } = ( %$self, @_ );
    return $self;
}

sub cp_stats {
    my $self = shift;
    my ($values, $ranges, $stats) = @_; # optional stats
    error("no values given") unless (ref($values) eq 'ARRAY' && @$values > 1);
    my ($distrib, $mean, $median, $stdev, $samples, $valid);
    ## total sample dies
    $samples = scalar @$values;
    ## compare mode
    my $comp_mode = defined($self->{Comp_Mode}) ? $self->{Comp_Mode} : 0;
    ## get distribution
    if ($comp_mode & 4) { # discrete values distribution
        $distrib = discrete_distribution($values, $ranges);
    } else {
        error("invalid ranges array") unless (ref($ranges) eq 'ARRAY' && @$ranges > 1);
        $distrib = numerical_distribution($values, $ranges, $comp_mode);
    }
    ## ASSERT: description items === distribution items
    error("number of items mismatch: desc: $#{$distrib->[0]}, count: $#{$distrib->[1]}")
        unless ($#{$distrib->[0]} == $#{$distrib->[1]});
    ## valid dies (those with valid data)
    if ($distrib->[0][ $#{$distrib->[0]} ] eq 'n/a') { # "n/a" items found
        $valid = $samples - $distrib->[1][ $#{$distrib->[0]} ];
    } else {
        $valid = $samples;
    }
    ## get statisticals
    if (defined($stats)) {
        for (qw/Median Mean Stdev/) {
            error("'$_' in \$stats structure not numerical: $stats->{$_}")
                unless (Statistics::is_numerical($_));
        }
        ($median, $mean, $stdev) = @{$stats}{qw/Median Mean Stdev/};
    } elsif ($comp_mode & 4) { # discrete values have no statistics
        $median = $mean = $stdev = 'n/a';
    } else {
        my $valid_values = ($valid == $samples) ? $values :
            [ grep(defined($_) && $_ ne '', @$values) ];
        #message("valid values is: $valid_values (count: @{[scalar(@$valid_values)]}/@{[scalar(@$values)]})");
        $median = Statistics::median($valid_values);
        $mean   = Statistics::mean($valid_values);
        $stdev  = ($valid == $samples) ?               # whether all sample dies' data are valid
            Statistics::stdev($valid_values, $mean) :  # the entire population
            Statistics::stdev_s($valid_values, $mean); # a sample of the population
        for ($median, $mean, $stdev) { $_ = 'n/a' unless (defined($_)); }
    }
    return +{
        Median       => $median,
        Mean         => $mean,
        Stdev        => $stdev,
        Sample_Dies  => $samples,
        Valid_Dies   => $valid,
        Distribution => $distrib,
    };
}

### class method
sub numerical_distribution {
    my $argc = grep(defined, @_);
    my ($values, $ranges, $comp_mode) = @_;
    error("only accepts 3 arguments") unless $argc == 3;
    ## compare subroutines
    my ($comp0, $comp1) = @{$comp_subs[$comp_mode]};
    my @distrib;
    ## sort to different ranges/categories
    my ($extra, $na) = (0, 0);
    my @counts = map 0, 0 .. $#$ranges;
    ## COMPARE DESCRIPTION
    ## sample range: [1, 2, 3, ]
    ## type-1: mode: 2'b00 or 2'b10
    ##         'greater-than' mode: { <1 }, >=1, >=2, >=3
    ##         COMPARE: do { <1 } first, and then traverse BACKWARDS!
    ## type-2: mode: 2'b01 or 2'b11
    ##         'less-than' mode: <=1, <=2, <=3, { >3 }
    ##         COMPARE do { >3 } first, and then traverse forwards
    ## COMPARISON STEPS
    ## 0. according to the "compare mode", decide the following things:
    ##    a. decide range edge, whether left (0) or right ($#$ranges)
    ##    b. decide comparing order when doing range sorting, i.e.,
    ##       from left to right (0 .. $#$ranges) or reversed
    ## 1. if $val is not defined, or empty, then it's deemed as "na" (N/A)
    ## 2. if $val is on the edge (eg, <LEFT or >RIGHT), then it's "extra"
    ## 3. else, sort to specific range
    ## 4. otherwise, error
    my $index_extra = ( $comp_mode & 1) ? $#$ranges : 0;
    my @indices = 0 .. $#$ranges;
    @indices = reverse @indices if !( $comp_mode & 1);
each_value:
    for my $val (@$values) {
        if ( !defined($val) || $val eq '' ) {
            $na += 1;
        } elsif ( $comp1->($val, $ranges->[$index_extra]) ) {
            $extra += 1;
        } else {
            for my $i ( @indices ) {
                if ($comp0->($val, $ranges->[$i])) {
                    $counts[$i] += 1;
                    next each_value;
                }
            }
            # still not sorted!
            error("value: $val: range: [@$ranges] at mode '$comp_mode': not sorted!");
        }
    }
    @distrib = ($comp_mode & 1) ? (@counts, $extra) : ($extra, @counts);
    push @distrib, $na if $na;
    my $desc = CPStats::limit_desc($ranges, $comp_mode, $na);
    return [ $desc, \@distrib ];
}

### class method
sub discrete_distribution {
    my $argc = grep(defined, @_);
    error("only accepts 1 or 2 arguments") unless ($argc == 1 || $argc == 2);
    my ($values, $ranges) = @_;
    error("seems there are some numerical values?")
        if ( Statistics::is_numerical( $values->[0] )
             || Statistics::is_numerical( $values->[$#$values] ) );
    my ($na, %distrib);
    for (@$values) {
        if (!defined($_) || $_ eq '') {
            $na += 1;
        } else {
            $distrib{$_} += 1;
        }
    }
    my (@desc, @distrib_ret);
    # if the user's specification of possible values
    # does not include all values given, we'll signal
    # it in an extreme way
    if (defined($ranges)) {
        my %given_vals = map {; $_ => 1;} @$ranges;
        for (keys %distrib) {
            error("found: $_: but not present in given ranges: [@$ranges]")
                unless ($given_vals{$_});
        }
        @desc = @$ranges;
    } else {
        @desc = reverse( sort(keys %distrib) );
    }
    @distrib_ret = map($distrib{$_} ? $distrib{$_} : 0, @desc);
    for (@desc) { # try to use a more friendly description string
        if ( defined( $discrete_desc{$_} ) ) {
            $_ = "$discrete_desc{$_} ($_)";
        }
    }
    if ($na) {
        push @distrib_ret, $na;
        push @desc, 'n/a';
    }
    return [ \@desc, \@distrib_ret ];
}

1;

__END__

=head1 NAME

CPStats - CP Statistics calculation module

=head1 SYNOPSIS

C<my $cps = new CPStats::( Data =E<gt> $wafer_data );>

C<my $stats = $cps-E<gt>calc_stats();>

C<print "Hello, world\n";>

C<my %hash = Name =E<gt> Xia Zuoling, Gender =E<gt> 'male';>

For your reference, before you start to write codes, you may
want to have a look at the section
L</Definition of Standard Deviation> first.

=head1 DESCRIPTION

=head2 Definition of Standard Deviation

The following content is from B<Wikipedia (R) - the free encyclopedia>

The link is:

I<L<http://en.wikipedia.org/wiki/Standard_deviation>>

And since our great government banned
L<http://wwww.wikipedia.org>, here is the alternative links:

=over 1

=item *
I<L<http://encyclopedia.thefreedictionary.com/standard%20deviation>>

=item *
I<L<http://www.answers.com/standard%20deviation>>

=back

The standard deviation of a random variable  X is defined as:

             _______________________         _______________________
            /   /             2  \          /                    2
sigma  =   /  E(  ( X - E(X) )^   )   =    /  E( X^2 ) - ( E(X) )^
         \/     \                /       \/

where E(X) is the expected value of X.

For discrete random variable X:

E(X) = mean(X) ==> X'

                     __N__
                 1   \             X1 + X2 + ... + Xn
X' = mean(X) =  ---   >    Xi  =  ---------------------
                 N   /____                 N
                      i=1

So the standard deviation of a discrete uniform random variable X is
calculated as follows:

=over

=item 1
For each value Xi, calculate the difference (Xi - X') between Xi
and the average value X'.

=item 2
Calculate the squres of these differences.

=item 3
Find the average of the squared differences. This quantity is the
variance (sigma^2).

=item 4
Take the square root of the variance

=back

In real world, for a sample of values X1, ..., Xn from some larger
population, a "sample (or estimated) standard deviation" is defined as:

          _____________________________
         /        __N__
        /    1    \      /         \ 2
s  =   /  -------  >    (  Xi - X'  )^
     \/    N - 1  /____  \         /
                   i=1

The reason for this definition (N-1 instead of N in the denominator)
is that I<s> is not an "unbiased deviation", unless the population is
normally distributed. This definition of standard deviation tends to
overestimate a little, i.e., this calculated "stdev" value is a little
larger.

And the following definition, contrarily, tends to underestimate the
real sigma:

          _____________________________
         /        __N__
        /    1    \      /         \ 2
s  =   /  -------  >    (  Xi - X'  )^
     \/      N    /____  \         /
                   i=1

This form has a uniformly smaller "mean square error" than does the
"unbiased estimator", and is the "maximum-likelihood estimate" when
the population (or the random variable X) is normally distributed.

=head1 BUGS

A lot of bugs exist!

=head1 AUTHOR

=over

=item *
Xia Zuoling <I<xiazuoling@163.com>>

=back

=head1 COPYRIGHT

=over

=item * `CPStats' I<ver. 1> Copyright (C) 2007

=item * POD update: I<2007-4-30 15:39:12>

=back

=cut

# vim: set filetype=perl:
